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ABSTRACT 


The generation and propagation of small amplitude 
acoustic waves in a homogeneous, loss-free, compressible 
fluid is studied by the finite element method. A diaphragm 
mounted in an infinite rigid baffle generates acoustic 
waves in a semi-infinite fluid region. Steady-state pres- 
Pee distribution as found for a hemispherical region with 
boundary reflection suppressed through use of a radiation 
condition. The computer program developed for the purpose 
utilizes iso-parametric finite elements with curvilinear 
boundaries. incorporated in the program is a versatile 
Hieom £6eneracor whieh minimizes the quantity of input data. 
Macepuvaple agreement with analytic results is obtained when 


there are at least four elements per wave-length. 
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i. INPRODUCTION 


Tie “fLpeld lok acousities plays an ampertvans Grole in the 
development of the complex sonar systems used in the United 
States Navy today. It is theLaincent oof Chislparper tie, open 
a new avenue to the study of the active sonar system. 
Beuvenvtion is confined to. the transducer and the: periodic 
pressure field induced in a fluid medium by harmonic 
Sserllations of the transducer piston or diaphragm. 

the study as restricted to the steady-state performance 
of a transducer mounted in an infinite rigid baffle and 
Fadiating into a half-space filled with a homogeneous fluid. 
Both the pressure field and the acoustic impedance of the 
Gransducer are determined. Figure 1 illustrates the 


geometry of the region. 


Fluid Medium 


Diaphragm 


Figure 1. Geometry of Region. 


The method used in this study is the two-dimensional 
finite element technique. Since it is manifestly impossible 
GOsrepresentcan inimite rection Wit acelin be snumber iof 
finite elements, the stratagem of the anechoic chamber is 
Simulated mathematically. This is accomplished by using 
radiation (non-reflection) boundary conditions. 

the diaphragm, being circular, forms a surtace of 
revolution when deflected. Thus, the instantaneous pres-— 
Sure distribution has rotational symmetry about an axis 
teeough the diaphragm center and normal to the rigid bafilte 
Geeexis of Fisurel), <A finite element assembly representing 
mee refion Of anberest may be constructed by dividing the 
eorresponding portion of the r-Z plane of Figure 1 into sub- 


areas as shown in Figure ed. 


r Non-reflection Boundary 
/ (indicated by dotted line) 
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Diaphragm 
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Figure 2. Finite Element Grid 
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In Figure 2, an individual rectangular element represents 
bhe cross-section of a toroidal element @enéeraced by revolu— 
Giron about the, Z axis;  ThemotationaY symmetry permics 
calculation of element properties using two-dimensional 
technique by simply taking the thickness (normal to the 


plane of the figure) proportional to r. 


IT. FINITE ELEMENT FORMULATION OF 


THE GOVERNING EQUATION 


Pe FQUATION OF FLUID MOTION 
The governing partial differential equation for the pres- 
sure p in a homogeneous, inviscid, compressible fluid as 


sed in acoustics (e.g., see Ref. [1]*) is 


Vp =P Ge 
C 


where 7 is the haplacian @peratdr, ¢ is the acoustic 
moc ty, and Superior dots represent diiferentiation with 
maopect, tO times, For ‘the purpose,.of the present .study, 


boundary conditions take the form 


ap 


The dees = 


where the coordinate n is measured along the normal to the 
boundary surface, p is the fluid density, and Ne is the 


Pemponent of fluid particle velocity in the n direction. 


B. SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS 

The finite element method replaces Eqs. (1) and (2) by 
a set of simultaneous ordinary differential equations with 
independent variable time and dependent variables consisting 
of pressures at a large number of discrete points (nodes) 


of the region. The region is divided into a number of 


*\umbers in brackets designate references listed on page 68. 


sub-regions called "elements." Within each element it is 
postulated that the pressure p may be expressed in terms of 
the pressure Ds at the n nodes of the element. This rela- 


tion is expressed as 
n 
ee PNA oe (3) 


where the "shape" (or interpolation) functions N, are 
funeGions of the spatial coordinates. If the nodes are 
foeated at suitably chosen points on the element surface, 
unique sets of shape functions will guarantee the necessary 
continuity of pressure p along interelement boundaries as 
pounted out: by Zienkiewiez [2]... The elements used in this 
study are two-dimensional eight-noded iso-parametric 
elements. Further details are given in Appendix C. 

It has been shown by Zienkiewicz and Newton [3] that 
the finite element equation for k nodes, m of which are 


Heocated on the fluid-diaphragm interface, may be written as 
as : ik te 
Olin + [Dito + [alipt ==. pbb) 16) (4) 


where {p} is a k xl vector of nodal pressures, [Q], [D], and 
baleare k xk symmetric matrices, {6} is an mx I vector of 
diaphragm displacements, and seat 1Sa xm matrix. Elements 
of [0], [D], [H], and [L]_? are all real constants. Formulas 
for calculating these matrices are given in Appendix E. The 
terms [D]{p} and ll te) represent boundary conditions and 


moet! be discussed in following sections. 


C. DIAPHRAGM BOUNDARY CONDITION 

The diaphragm displacement 6 is a prescribed function of 
time. It can be represented by choosing a number of points 
(nodes) on the diaphragm and using a set of shape functions > 


N, to interpolate the local displacement as 
é m 
& Saas? (5) 


where the - are the displacements at the m diaphragm nodes. 

The boundary condition of Eq. (2) is formulated in terms 
of the normal component ee Ont 1 UC" par pase ihe’ vedo Than 
the diaphragm surface this must, in the absence of cavita- 
tion, be equal to 6. Using this equality, the fluid- 


diaphragm interface—condition gives the term 

-p{L] {8} 
in the finite element equation (Eq. (4)). The resulting 
Wector has nonzero elements only for those pressure nodes 
iecaved on the diaphragm surface. it is noted that the 


peundary condition on the rigid baffle is a special case, 


in that the particle velocity v, = 0, which implies a = 0. 


D. RADIATION (NON-REFLECTION) BOUNDARY CONDITION 

A normally incident wave meeting the interface between 
two media will not be reflected if the specific acoustic 
impedances of the two media are equal. The specific acoustic 


impedance z is defined to be the ratio of pressure to 


10 


perticle velocity. Hor plane harmonic. waves, unts as a. meeaal 


Quantity and the pressure can be given as 


Do ravlve ie (6) 


for a wave progressing in the positive m direction. Differ— 
entiating Eq. (6) with respect to time and using this result 


to eliminate ae from the boundary .conditien Giquiiee ) jmeawes 


opps 


a . 
on @ 


p CO) 
Equation (7) then becomes the boundary condition for non- 
reflection of normally incident plane waves. The term 

[D} {p} 


mong. (4) results from this criterion. This vector has 
nonzero elements only for those nodes which lie on the non- 


fersecting portion of the boundary. 


aait 


LUT. (‘GEOMETRY OF TRE REGION 


The initial solution of the problem using the geometric 
Sentifuration of Figure? resmiiied in) sient rveant error an 
the pressure distribution when the solution was checked 
eaeaimst known results. The principal cause of this error 
mas believed to be the shape of the non=-reflection boundary. 
By the nature of the geometry, the harmonic oscillating 
Giaphragm will create progressive pressure waves which | 
propagate outward and strike the boundary of Figure 2 at 
Various angles of incidence. The particle velocity, Vis 
Will not always be normal to the non-reflection boundary 
and the boundary condition required by Eq. (2) will not be 
Baersried. Thus), a modification to the representation of 
tae boundary to provide substantially normal incidence led 
momene feometry illustrated inFigure 3. For this modified 
Paciation boundary, in three dimensions the surface is 
hemispherical. 

in order to represent adequately the circular arc of 
the non-reflection boundary, two-dimensional iso-parametric. 
finite elements are used. The element thickness normal to 
meme plane of the figure is directly proportional. to the 
radius r. Iso-parametric elements are further discussed 
mm Appendix C, 

Note the distortion Of the grid in Figure 3. The reason 


for this distortion and the method of generating the grid 


L2 


x 
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Diaphragm | 


Figure 3. Modified Finite Flement Region 


pattern is discussed in Appendix D. The distorted 
"rectangles" still represent the cross sections of toroidal 


elements. 


1S, 


IV. WAVE THEORY AT RADIATION BOUNDARY 


A one-dimensional study of a pulsating sphere acting on 
eeriuad medium resulted in @ more aceurate solution when 
Spherical wave theory at the non-reflection boundary was 
maeroduced. Based upon the success of these results, a 
further modification was made to the computer program 
(Appendix B) to incorporate spherical wave theory in the 
two-dimensional study. This was accomplished by developing 
the non-reflection boundary condition for spherically diverg- 
ing waves as shown below. 

The specifie acoustic impedance for spherically diverg- 


ing waves is given in Ref. [1] as 
Z = a= = 0G === == (8) 
rm : 


Mere wW is the circular frequency and r is the radial 
distance from the "center." The particle acceleration Vn 


can be written as 

as SVE, (9) 
which allows the boundary condition to be represented by 
(a ree 

or PJOVA, (10) 


From Eqs. (8) and (10), the boundary condition is 
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The imaginary term on the rine hand Sige On. Pa. wile)aats 


mepresented by the term {[D]{p}.= jw[DIl{p} in Eq, (4). 


the contribution of the real term can be written as 


c 
pel {p} 


Maybe 


(Cs) 


maga can be added to the system of equations given in Eq. (4%. 


The spherical contribution is only added to the nodes 


which lie on the non-reflection boundary, thus providing a 


spherical wave approximation to the plane wave theory at 


Gae non-reflection boundary. 


iD 
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V. COMPUTER SOLUTION TECHNIQUE 


Assume the diaphragm is oscillating with a harmonic dis- 


placement represented in the complex form by 


eh eae 7 (1) 


where HOR is a vector of prescribed displacement amplitudes 
at the m nodes of the diaphragm-fluid interface. Elements . 
of Ee) are real, | The presume in the fluid medium aseemen 


represented by 
{p} = {p yet" (sy) 


where ip is a kxi1 vector of pressure amplitudes (complex) 
at the k nodes of the region. 


inservtaneatgse. (14). andy (15). into Eq. s( 4) eives 


([H] - w°[Q] + jwID]){p} = w°{B} (16) 


where {B} = a oe). 

In the computer program, the system matrices [H], [Q], 
[D], and {B} are assembled from the contributions of each 
migividual findte element in various subroutines. Since 
only a small number of nodes are located on the non- 
ferlection boundary, the matrix [D] is sparse. In order to 
reduce computer storage requirements, [D], which is symmetric, 
@™= Scored as three vectors: the principal diagonal and the 
two diagonals immediately above. At the same time, the 


approximated spherical contribution 10] at the 


16 


non-reflection boundary is computed. These quantities are 
real and are added directly to the wees nn in the appro- 
prrave Locations. 

The final subroutine introduces the value of omega and 
multiplies the matrices [Q], [DJ], and {B} by the appropriate 
power of omega before the left-hand side is summed to form 
one complex matrix. Thus the system of equations is 


reduced to 
[Sitohe=oikt | Cia) 


meere [5] 1s ak xk complex matrix and {Ff} is a column 
maecor Of Order ks | in whiten only m elements are nonzero. 

The complex pressure amplitudes in Eq. (17) are found 
by using the Gaussian elimination and back substitution 
method. 

The absolute value of the pressure and the phase angle 
pee caliculated ian the program and printed ous in array 
form. Additionally, a system energy balance and the 
diaphragm acoustic impedance are determined. 

The program provides for nappen pee. cet eh solution by 
using incremented values of omega so that a range of 


Frequencies can be studied. 


Lif 


Vi. SRESULUS: 


Problems studied were.chosen tor vers iy, vmae accuracy .or 
the method through comparison with available, exact 
(analytic) solutions. Results for each problem are 


Giscussed below. 


Beeuvemn i Srisvon in a higidy lube 

To verify the correctness of the computer program, 
finite element solutions were obtained for the problem of 
aearmont eallysoscilllating eireular pisvom at’ one end of 
a fluid-filled, rigid pipe of equal diameter. ‘An infinite 
length was simulated by terminating the finite element 
region with a plane non-reflection boundary at a distance 
Perene piston diameter. This is, of course, a one dimen— 
sional problem with uniform pressure amplitude throughout 
and phase angle is a linear function of distance along 
the pipe. 

Detailed results are not presented for this essentially 
maovral test problem. “it was found, however, that there 
was excellent agreement with theory for long wave-lengths. 
For wave-lengths less than approximately four times the 
element width (measured parallel to the direction of wave 
propagation), the accuracy deteriorated. This is under- 
standable, since the shape functions employed represent 
the instantaneous pressure profile across an element by a 


parabolic are. The lower limit on wave-length corresponds, 
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for a givenvelement gerade tor an upperilLini tom irequency - 
Eeob Lem, 20. Pustonvinvan, Tai inv emsann ie 

The problem explored most extensively is that of a 
harmonically oscillating circular piston surrounded by an 
Pmianite raigidibaifle and radiating antova wema—-intinite 
Pigidjregion. Yor this problem, thervradiavdon beundary is 
a hemispherical surface with the center located at the 
piston center and the radius is four times that of the 
mee tbon. The element grid, (see'\ Figure: 4, Appendix’ A) contains 
16 elements. It is derived by mapping a 4 x4 square array 
imome Che = .joplane into a*quadrant of a-cirele. 

Analytic results for this problem are readily avail- 
moeuteie.) iseeatiivr [4] )oondn Table I, the pressure ampli- 
tudes and phase angles along the Z axis, obtained by the 
finite element method (FEM), are compared with exact 
results. Two sets of finite element results are given. 

For those designated FEM1, the spherical wave correction at 
eer radiation boundary is included. For FPEM2, the correc-— 
Peon 1s omitted. Details are given in Table I for two 
circular frequencies, w = .1 and .6 rad/sec. The corres- 
ponding wave-lengths are approximately 63 ft and 10 feet. 

Examination of Table I shows that FEM1 results give 
generally good agreement with exact results. Those for 
FEM2 are considerably poorer, although the Gs soecueneiee 
are smaller for w = .6 rad/sec. This is understandable since 
the i7jrvecentribution to specifresacoustie |\impedance is 


relatively less important as w increases (see Eq. (12)). 


a) 


It is believed that a more refined finite element mesh 
(more elements) should afford significant improvement. 
Since the computer program for this problem required 
320K bytes of core storage, it was not feasible to 
MmvesLigare thas without major program revisions. 

In predicting, transducer, performance ,uhe, resultant fituad 
reaction on the piston is needed. This is normally reported 
ai the form of ~ piston radiation impedance, which is- complex 
ena oS found by dividing the piston foree by the piston 
velocity. The real-part of the radiation impedance is 
called the radiation resistance and the imaginary part is 
Gered the radiation reactance. 

Table IL provides a comparison with exact values for 
the radiation resistance and reactance. Results are given 
ieee! and BEM over a range of frequencies, Again it is 
Opeerved that the spherical correction significantly improves 
the results. Note that the radiation resistance obtained 
by PEML continues to give good results through © = .9 rad/ 
See 

Computer solutions of this problem were obtained for a 
range of frequency from w = .1 rad/sec to w = 1.2 rad/sec 
by .1 rad/sec increments. Comparisonsat the upper end of 
this range are omitted from Tables I and II because the 
Geparture from exact results does not permit useful employ- 
ment of the finite element method without grid refinement. 


Complete computer results for pressure amplitude, phase 
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TABLE I 
PRESSURE AMPLITUDE AND PHASE ANGLE ON AXIS 


FEML FEM2 FEM] FEM2 


ft 1p /Ge inte deg 
= 0 02061 | .01794 7. \e eaten 
i 2 +00889 | .00676 Ml 23ihe 
rat \ 00469 | .00376 eenGee 
ny 6 00329 | .00337 beac om 
4 8 “00247 | .00334 ol eae s 
56 0 625 eas. 
6 2 DB ite| =86..2 
.6 y Suilten 6 =a SOS 
% 6 106 Lo | 22a6eG 
6 8 .084 <2 | =2090R4 
TABLE IT 


PISTON RADIATION RESISTANCE AND REACTANCE 


SS SSE SSS <i = = EE — 


REACTANCE 


EXACT FEMI FEM2 


FEM1 FEM2 
_ ibesec | 1b .sec | wiby sec Lbesec lb see 
ft ft 


eel 249 e254 2h ae 1.68 
eee 29 .978 h, Bee Sig 56 
#5 eae eueke Die Ds boos 
4 Siaiee By O2 (ee Gran cleus 
25 Bee 5.34 oe {. 7.43 
a6) (One (lee Oe he 00 
af 8.89 uate oe 8. Ohi 
a8 10S5 1 LOS 53 On (as 02 59 
9 n.90 TAR NF (e ie 6.69 


FEM1 - Finite element results with spherical correction. 
FEM2 - Finite element results without spherical correction. 


Phase angle given is for pressure relative to piston 
displacement. 


Fiuid density = 1 ib cases tt: Pistom Radius = 2 ft. 
Peoustic Velocity = 1 ft/sec. Piston Jmpdatude = 1 wt. 


eae 


angle, and impedance components at several frequencies are 
Ziven in Appendix A. 
imic@be Ciiy=3 40h arabodieeDiLephraam ranma neem sa iiS 

This problem, ‘a modification of -the-previious—probi-em, 
meplaces the piston by a circular diaphrarm having @ faxed 
edge and having a deflection surface in the form of a 
paraboloid of revolution. Study of this problem was under- 
taken in an effort to improve agreement in the energy 
balance check (discussed below) incorporated inte the 
peoeram. It was believed that elimination of fluid 
melocity discontinuity.at the edge of the piston in problem 
fmrecht resuie an a significant improvement. 

Analytie results for this problem are not available in 
Svencgard acoustic texts, but conventional techniques can 
be utilized to obtain the pressure p along the Z axis in 


complex form as 


i, 
a3 ae -j==-] , 
p= Pipe So fis 2 eb =, (eee aa Bae Mee Cre) 
Wa (AS 
where 
Sg = displacement amplitude at diaphragm center 
a = diaphragm diameter 
i "as E 2 2% 
Ss = “sien distance — ©Z -¥ a) 


For this problem only a single finite element solution 
fyEM3) was obtained. The spherical correction is included. 


Mable Ill gives a comparison with exact results (Eq. (18)) 


oie 


ils ah > Layer 
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a = _ ne = * 
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PRESSURE AMPLITUDE AND PHASE ANGLE ON AXIS 


| pRessuRE AMPLITUDE PHASE ANGLE | 
w z EXACT pem3 | EXACT FEM3 
rad/sec ft 1b/ft? lb/ft? deg 
a 0 Bg 322 01127 Sah 
Gi 2 00438 00461 212u6 
a i 00240 002H2 EA 
eo. 6 .00164 .00166 — 35.5 
an g 00124 00125 hi 
Sao 6 4598 4934 23.8 
6 2 565 1590 S755 
6 \ .0864 YOS'6 7 ~142.3 
6 6 .0589 .0584 -209.4 
6 8 0445 0436 -278.1 


FEM3 - Finite element. results with spherical correction. 
Phase angle is for pressure relative to diaphragm displacement. 
Pituid density = 1 Mp eeeyre 

Meoustic velocity = 1 ft/sec. 

Daaphragm radius,.= 2 ft. 

Diaphragm amplitude (center) = 1 ft. 


for pressure amplitudes and phase angles for values of 
® = .1 and .6 rad/sec. 

Examination of Table III shows that agreement near the 
faaphragm (small values of 2) is not as good as FEMI in 
Table I. At the radiation boundary however, the FEM3 errors 


are comparable to those for FEMI. 


Energy Balance Check 


Prior experience, by others, suggested that an energy 
balance check would be a useful feature of a finite element 
program. Specifically, since solutions represent steady— 
state behavior, the input power, averaged over a cycle, 
sheuld equal the corresponding average value of the output 
pewer, Input occurs* at the piston, or diaphragm, and output 
at the radiation boundary. Expressions for these quantities 
ere derived and incorporated adn the computer program. 

A number of program errors in the energy balance calcu- 
lations (Appendix B) prevented valid results until the 
finale stagesorethe investigation:: Thus, the values of 
average power in and average power out are not included in 
the computer output (Appendix A). Recent computations, 
fae propram errors corrected, have shown excellent: agree— 
ment of these quantities to at least six significant digits 
weeune lower range of the circular frequency. Even at 
higher frequencies, where solution accuracy (pressure 
amplitudes and phase angles) Geteriorates markedly, five 


Significant digit agreement is maintained. 


eu 


VII. CONCLUSIONS 


It has been shown that. the finite element method ican 
provide acceptable results concerning generation and prop- 
agation of acoustic waves. The parabolic iso-parametric 
Slenent provides .an economical, accurate means of repre= 
senting regions having curved boundaries. 

The principal’ Limitation ef the method is evident. from 
the empirically deduced results that, for acceptable 
accuracy, a single element must not span more than one- 
quarter wave-length. For wave-lengths of the order of 
magnitude of the piston radius, the number of nodes 
required would be around 900. For the present program, 
which already utilizes 320K bytes of core storage for 65 


Mee@es. this is far beyond capabibity. 
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APPENDIX A 


COMPUTER OUTPUT 


The rows of _ pressures and phase angles shown in the 
computer output correspond to the numbering system shown on 
the finite element region, Figure 4. 

The zeros printed in the array of pressures and phase 
angles do not indicate computed values, but facilitate 


Drtnoing the output in array form. 


Row 1A 
Row 20 || 
Row 3 M20 
Row 4/28 29 


Row 37 38 


Row 
Row 
Row 


Row 


Figure 4. Format for-Pressure and Phase Angle Array. 
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The output data are arranged according to the following 


sequence. 


REMI, Piston in infinite batile with spherical correcunen 
added at the non-reflection boundary. 

reve. Piston imeinfinite baffle without sphemical .cormec— 
tion added at the non-reflection boundary. 

FEM3. Diaphragm in infinite baffle with spherical correc-— 


tion added at the non-reflection boundary. 
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APPENDIX C 


ITSO-PARAMETRIC ELEMENTS 


The family of iso-parametric elements, recently developed 
by Zienkiewicz, Irons, and others, has brought new versa- 
tility and power to the finite element method. The special 
usefulness of these elements is that they may be shaped to 
eoniform vo curvilinear boundaries. Attention is confined 
here to two-dimensional parabolic elements. A full descrip-— 
muon Of the enbhire family is e@iven in Ref. [2]. 

Consider vene rectvangulbar element of Higuce 5. The 
element has four corner nodes and four midside nodes, 
membered as Shown. Points in the interior, or on the 
boundary, may be located in the local (&,n) coordinate 
system. The dependent variable, pressure p in the present 
geplication, is defined at any point by Eq. (3), repeated 


here 
N.DP, C32 


where 
N, = (lt nn,)(1 + €6,) (nn, + €&, -1) 


bor ‘corner node 1, and 


MS s(1- ete ~ ee! alll tees) 
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Paeure 5. Element Configuration in Local’ &.n Coordinates 


for mid-side node 1, 


(E, on; ) are coordinates of node i, 


and Ds is the pressure at node i. 
Note that each shape function Ny has the value unity at node 
i and is Zero at all other nodes. 
A curviilimeanr- element, such as CpREnOE Figure 6, may 
be produced by mapping the rectangle of Figure 5 into the 


may Dulane: by the relations 
8 
Spel NO dal Sy Ci) 


where (X, 594) are the cartesian coordinates of node i. 
The shape of the element in the x,y plane is uniquely 


determined by specifying the coordinates (x, 54) for the 
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1 


Figure 6. Curvilinear Elements Mapped 
in (Cartesian, Coordi naves: 


~ 


eieht nodes." The edges of this curvilinear element are 
parabolic. 

Adjacent elements "share" the same three nodes at their 
interface. When mapped to cartesian coordinates, these two 
adjacent elements then have the same curvilinear shape and 
the same x and y coordinates at nodes. Thos provides con— 
tinuity in the pressure distribution across element 
boundaries. | 

The region developed for the present study necessitates 
the use of iso-parametric elements to provide the ei rcular 
non-reflection boundary. This boundary is constructed 
using parabolic are segments which are mapped by use of the 
shape Frunevions to "fit" a quarter circle boundary in car- 


4esian coordinates. 
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APPENDIX D 


FINITE ELEMENT MESH GENERATOR 


two subrouwGiness | aneluded: tn the tconputersprosram von 
Appendix B, should prove useful inv variety tof finite 
element problems. He EENE ty these subroutines subdivide 
the region, systematically assign element numbers and 
hnede numbers, calculate nodal coordinates, and produce 
a computer plot of the region to be studied. As an example, 
Consider the megion of Figure 7. 

The region is taken to be one large iso-parametric 
element (a mapping of thie rectangle of Figure 5). Its shape 


is prescribed by specifying the (x,y) coordinates of the 


Figure 7. Basic Geometry of Finite Element Region. 
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Vee is)? 


eight nodes. Subdivision of the element can then be 
accomplished at predetermined values of —€ and 7 by using 
Equation (Cig )e 

Subroutine NODE (Appendix B) requires two input 
patcemeters , the number of divisdions Gesiredvamyehe x diaee= 
mon and the number of divisions desired in the y @irect ion, 
These parameters are eee to calculate the increment values 
of € and n for subdividing the element into smaller 
elements. The subroutine starts numbering global nodes at 
the upper left hand corner (&=-1, n=1) then adds Amevement 
values of € and assigns global node numbers from left to 
meeht on the boundary n=l. Similarly, arpeer, imcrementing 
MH, the second row of global node numbers is assigned, etc. 
Global nodé numbers are assigned to all corner nodes, mid-= 
side nodes and the center of each subdivided element (9 
nodes per element). The center node provides unique 
aoenvification for the element and simplifies the printing 
of computer output in a readable form. Subsequently all 
modes are renumbered, omitting the center nodes. 

The cartesian coordinates corresponding to the global 
node numbers are calculated by a similar process in the sub- 
routine XYFORM and are used in the program to plot the 
finite element grid. Figure 8.illustrates a computer plot 
of the finite element mesh generated for a quarter circular 
region (Figure 7) with an eight foot radius. The number 
Ot subdivisions im both the x end y direction Ws equal ve 


four. The plotting routine produces straight segments 


Figure 8. Finite Element Mesh Generated 
by Subroutine XYFORM. 


between nodes. In this instance, the element boundaries 
gee all parabolic arcs. except, for those dying om the x and 


y axes. 
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APPENDIX E 


MATRIX ELEMENT FORMULAS 


individuals estements Of; the matriuecs 10 li. (oieeslh hemmed 
tele ot the system of ordinary differential equations, (Eq. (4) 


ere given by Ref. [3] as 


i 

di: » See Se Ne Nd 

ie es 

Ae ee Nas 

ay} e Se aes 
{= ihn Nn dS 

aL aes 

oN, 9N, 23N, aN, 3N, aN, 

hh ea ax By co emer 


where Re and Se are the element region and external boundary 
respectively. 

The shape functions Ny; Ne are specified in terms of 
fecal coordinaves and it is advantageous to perform the 
integrations in these coordinates. Gaussian quadrature is 


employed. Details of the procedure are given in Ref. [2]. 
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